function cT_k = fft_test
%%%%%%%%%%%%%about%%%%%%%%%%%%%%%%%%%%
% option pricing via fft
% please complete comments
% try to extend the code and include
% other models: jump diffusion, sv, etc.
%parameters
S0 = 100;% initial underlying price
K = 80;% strike
r = 0.05;% interest rate
q = 0.01;% divident rate
T = 1.0;% maturity
n = 12;% montly frequency
sigma = 0.2;% volatility
alpha = 1.5;% add comment
eta = 0.25; % add comment
%% option pricing via fft
[~, cT_km]=genericFFT(S0, K, r, q, sigma, T, alpha, eta, n);
% output: option price in stage 1
cT_k = cT_km(1);
end

function phi = generic_CF(u, S0,sigma, r, q, T)
%%%  Computes the characteristic function for BS model.  
mu = log(S0) + (r-q-sigma^2/2)*T;
a = sigma*sqrt(T);
phi = exp(1i*mu*u-(a*u).^2./2);
end

function [km, cT_km]=genericFFT(S0, K, r, q, sigma,T, alpha, eta, n)
%%% Option pricing using FFT (model-free)%%%
N = 2^n;
%step-size in log strike space
lda = (2 * pi / N) / eta;
% choice of beta
%beta = np.log(S0)-N*lda/2 # the log strike we want is in the middle of the array
beta = log(K); % the log strike we want is the first element of the array
%forming vector x and strikes km for m=1,...,N
%discount factor
df = exp(-r*T);
nuJ = linspace(0,N,N) * eta;
psi_nuJ = generic_CF(nuJ - (alpha + 1) * 1i, S0,sigma, r, q, T)./((alpha + 1i*nuJ).*(alpha+1+1i*nuJ));
km = beta + lda * linspace(0,N,N);
w = eta * ones(1,N);
w(1) = eta / 2;
xX = exp(-1i * beta * nuJ).*df.*psi_nuJ.*w;
yY = fft(xX);
multiplier = exp(-alpha * km)/pi;
cT_km = multiplier .* real(yY);
end